BIOSTAT_1

PM 2.5 Concentrations

plotclr = brewer.pal(9, "PuBuGn")
class = classIntervals(Obs_value, 
                       9, 
                       style = "fixed", 
                       fixedBreaks = seq(min(Obs_value),
                       max(Obs_value), 
                       length = 9+1))
colcode <- findColours(class,plotclr)
plot(Longitude, Latitude, pch=20, col = colcode, main = "PM2.5 Measumrments on Select Sates")
US(xlim = range(Longitude),
   ylim = range(Latitude),
   add = T,
   lwd = 1,
   col = "gray")

map = GetMap.bbox(lonR = range(Longitude),
                  latR = range(Latitude), 
                  size = c(640,640), 
                  maptype = "hybrid")
PlotOnStaticMap(map)
convert_points = LatLon2XY.centered(map,
                                    Latitude,
                                    Longitude)
points(convert_points$newX,
       convert_points$newY, 
       col = colcode, 
       pch = 19)

x = Longitude
y = Latitude
model = lm(log(Obs_value) ~ x + y + State)
surf_est_surface = interp(x,
                          y,
                          model$fitted.values)
plotvar = surf_est_surface$z[!is.na(surf_est_surface$z)]
class = classIntervals(plotvar,
                       9,
                       style = "fixed",
                       fixedBreaks = seq(min(plotvar),
                       max(plotvar),
                       length = 9+1))
colcode = findColours(class, plotclr)
image.plot(surf_est_surface$x, 
           surf_est_surface$y, 
           surf_est_surface$z, 
           col = plotclr,
           breaks = seq(1.5,2.5, length = 10),
           zlim = c(1.5,2.5),
           xlab = "X", 
           ylab = "Y",
           main = "Estimated spatial trend")
US(xlim = range(Longitude),
   ylim = range(Latitude),
   add = T,
   lwd = 1,
   col = "red")

g = ggplot(data, aes(log(Obs_value))) + geom_density(aes(group = State, fill = State), alpha = 0.4) + ggtitle("Estimated Density of PM2.5 Concentration") + theme_economist()
ggplotly()

n = nrow(data)
dist_matrix = matrix(0 , nrow = n , ncol = n)
diff_matrix = matrix(0, nrow = n , ncol = n)
coordinates = matrix(cbind(Longitude , Latitude), nrow = n , ncol=2)
dist_matrix = rdist(coordinates, coordinates)
for(i in 1:n){
        for(j in 1:n){
                diff_matrix[i,j] <- (1/2)*(model$residuals[i]-model$residuals[j])^2 
        }
}
plot(as.numeric(dist_matrix), 
     as.numeric(diff_matrix), 
     type = "p", 
     pch = 20,
     col = "black",
     xlab = "Distance",
     ylab = "Squared differences",
     main = "Variogram Cloud Plot for PM2.5 Concentrations")

bins = (0.1)*((1:10)^1.8)
dt = data.frame(cbind(x, y, (model$residuals)))
empirical_variogram = variogram(model$residuals ~ 1 ,
                                locations = ~ x + y, 
                                dt,
                                boundaries = bins)
plot(empirical_variogram,
     col = "black",
     type = "p",
     pch = 20,
     main = "Empirical Semi-Variogram")

directional_empirical_variogram = variogram(model$residuals ~ 1,
                                            locations = ~ x + y, 
                                            dt,
                                            alpha = c(0,45,90,135),
                                            boundaries = bins)
plot(directional_empirical_variogram, main = "Directional Empirical Semi-Variogram")

exp_veriogram_model = fit.variogram(empirical_variogram,
                                    vgm(psill = 0.03, "Exp", 1, 0.01),
                                    fit.method = 2)
sph_veriogram_model = fit.variogram(empirical_variogram,
                                    vgm(psill = 0.03, "Sph", 1, 0.01),
                                    fit.method = 2)
gau_veriogram_model = fit.variogram(empirical_variogram,
                                    vgm(psill = 0.03, "Gau", 1, 0.01),
                                    fit.method = 2)
mat_veriogram_model = fit.variogram(empirical_variogram,
                                    vgm(psill = 0.03, "Mat", 1, 0.01, kappa = 1),
                                    fit.method = 2)
plot(empirical_variogram, exp_veriogram_model, main = "Exponential Model")

plot(empirical_variogram, sph_veriogram_model, main = "Spherical Model")

plot(empirical_variogram, gau_veriogram_model, main = "Gaussian Model")

plot(empirical_variogram, mat_veriogram_model, main = "Matern Model")